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To  understand  the  effect  of  chaos  in  the  classical  limit 
on  quantum  behaviors  has  been  a  field  of  interest  and 
active  pursuit  [1].  This  field,  named  quantum  chaos,  finds 
applications  in  many  fields  in  physics  such  as  condensed- 
matter  physics,  atomic  physics,  nuclear  physics,  optics, 
and  acoustics.  Previous  works  on  quantum  chaos  focused 
almost  exclusively  on  non-relativistic  quantum  systems. 
A  fundamental  question  is  whether  phenomena  in  non- 
relativistic  quantum  chaos  can  occur  in  relativistic  quan¬ 
tum  systems  described  by  the  Dirac  equation.  This  field  of 
relativistic  quantum  chaos  [2]  is  of  particular  importance 
because  of  the  relevance  of  the  Dirac  equation  to  graphene 
systems  [3]. 

Recently,  the  remarkable  phenomenon  of  chaos- 
regularized  quantum  tunneling  has  been  uncovered  [4], 
where  classical  chaos  can  suppress,  significantly,  the 
spread  in  the  tunneling  rate  commonly  seen  in  systems 
whose  classical  dynamics  are  regular.  For  example, 
consider  the  system  in  fig.  1,  which  consists  of  two 
symmetrical  cavities  connected  by  a  one-dimensional 
potential  barrier  along  the  line  of  symmetry.  When  the 
classical  dynamics  in  each  cavity  is  integrable,  for  suffi¬ 
ciently  large  energy  the  tunneling  rate  can  assume  many 

(a) E-mail:  xuan.ni@asu.edu 


Fig.  1:  (a)  A  closed  quantum  system  of  arbitrary  shape 
in  two  dimensions,  (b)  the  corresponding  opened  system, 
and  (c)  symmetric  double  well  used  in  resonant  tunneling 
computations. 

values  in  a  wide  interval.  Choosing  the  geometry  of  the 
cavity  such  that  the  classical  dynamics  become  chaotic 
can  greatly  enhance  and  regularize  quantum  tunneling. 
Heuristically,  this  can  be  understood,  as  follows.  When 
the  potential  barrier  is  infinite,  each  cavity  is  a  closed 
system  with  an  infinite  set  of  eigenenergies  and  eigen¬ 
states.  Many  eigenstates  are  concentrated  on  classical 
periodic  orbits,  forming  quantum  scars  [5].  For  classically 
integrable  cavity,  some  stable  or  marginally  stable  peri¬ 
odic  orbits  can  persist  when  the  potential  barrier  becomes 
finite  so  that  each  cavity  system  is  effectively  an  open 
quantum  system.  Many  surviving  eigenstates  correspond 
to  classical  periodic  orbits  whose  trajectories  do  not 
encounter  the  potential  barrier,  generating  extremely 
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low  tunneling  rate  even  when  the  energy  is  comparable 
with  or  larger  than  the  height  of  the  potential  barrier. 
The  eigenstates  corresponding  to  classical  orbits  that 
interact  with  the  potential  barrier,  however,  can  lead  to 
relatively  strong  tunneling.  In  a  small  energy  interval 
the  quantum  tunneling  rate  can  thus  spared  over  a  wide 
range.  However,  when  the  classical  dynamics  is  chaotic, 
isolated  orbits  that  do  not  interact  with  the  potential 
barrier  are  far  less  likely  and,  consequently,  the  states 
associated  with  low  tunneling  rates  disappear,  effectively 
suppressing  the  spread  in  the  tunneling  rate. 

In  this  letter,  we  address  the  question  of  whether  chaos 
can  regularize  tunneling  in  relativistic  quantum  systems. 
To  be  concrete,  we  study  the  motion  of  massless  Dirac 
fermions  in  the  setting  of  resonant  tunneling  to  facili¬ 
tate  comparison  with  the  non-relativistic  quantum  case. 
To  solve  the  Dirac  equation  in  a  confined  geometry  is 
extremely  challenging,  particularly  due  to  the  difficulty 
to  incorporate  zero-flux  boundary  conditions.  We  over¬ 
come  this  difficulty  by  developing  a  numerical  scheme 
based  on  constructing  a  physically  meaningful,  Hermitian 
Hamiltonian.  Our  extensive  computations  reveal  unequiv¬ 
ocally  the  existence  of  surviving  eigenstates  that  lead  to 
extremely  low  tunneling  rates.  As  for  the  non-relativistic 
quantum  case,  making  the  cavities  classically  chaotic  can 
greatly  regularize  the  quantum  tunneling  dynamics.  To 
explore  the  practical  implications,  we  consider  resonant 
tunneling  devices  made  entirely  of  graphene  [3] ,  and 
calculate  the  tunneling  rate  for  different  energy  values. 
We  obtain  qualitatively  similar  results  as  for  massless 
Dirac  fermions.  One  unique  feature  for  both  the  Dirac 
and  graphene  systems,  which  finds  no  counterpart  in 
non-relativistic  quantum  tunneling  devices,  is  the  high 
tunneling  rate  in  the  regime  where  the  particle  energy  is 
smaller  than  the  height  of  the  potential  barrier.  This  is 
a  manifestation  of  the  Klein-tunneling  phenomenon  [6-8] . 
To  explain  the  numerical  findings,  we  develop  a  theory 
based  on  the  concept  of  self-energies  and  the  complex 
energy  spectrum  of  the  non-Hermitian  Hamiltonian  for  the 
“open”  cavity. 

Consider  the  situation  where  a  relativistic  Dirac 
fermion  is  confined  within  a  two-dimensional  double-well 
system  in  which  two  symmetric  cavities  are  weakly 
coupled  by  an  electrical-potential  barrier  placed  in 
between.  The  system  is  governed  by  the  Dirac  equation 
=  H'ljjit),  where  the  general  form  of  Hamiltonian 
is  given  by  H  =  c(cx  -  p)  + /3mc2 ,  and  ^  is  a  two- 
component  Dirac  spinor.  Assuming  stationary  solution 
^(t)  exp(— iEt/ti),  we  obtain  the  steady-state  Dirac 

equation  Hi/j  —  E'lp.  In  two  dimensions,  ct  =  cr  =  (aXj  ay) 
and  /3  =  az  are  choices  satisfying  all  anticommuta¬ 
tion/commutation  relations  in  Dirac/Lorentz  algebra  [9]. 
There  are  two  major  challenges  in  numerically  solving 
the  Dirac  equation:  i)  a  proper  treatment  of  boundary 
conditions  and  ii)  an  efficient  and  physically  meaningful 
discretization  scheme.  Firstly,  for  boundary  conditions, 
we  can  either  replace  the  mc2az  term  by  a  potential 


U(r)az  in  the  Hamiltonian  and  let  U(r)  go  to  infinity 
outside  the  domain,  or  use  the  vanishing  current  condition 
j  •  n  =  0,  where  n  is  the  boundary  surface  normal.  We 
obtain  R e{el6(j)/x)  =  0,  where  (j)  and  x  are  the  components 
of  the  Dirac  spinor,  -0  =  (0,  x)T,  and  0  is  the  argument  of 
the  surface  normal  n.  It  was  demonstrated  [10]  that  the 
condition  is  x/0  =  ^exP(^)-  When  an  external  potential 
V  is  present,  such  as  the  barrier  in  our  tunneling  problem, 
E  is  replaced  by  E  —  V.  Secondly,  for  a  massless  fermion, 
we  can  write  the  Dirac  equation  as 

[v(cr  •p)  +  V]'ip  =  Eip,  (1) 

where  we  have  replaced  c  by  v  for  more  generalized  cases 
(e.g.,  in  graphene,  the  Fermi  velocity  is  +p^l06ms). 
Regarding  discretization,  the  difficulty  issue  is  fermion 
doubling  [11,12].  We  have  developed  and  validated  an  effi¬ 
cient  discretization  scheme,  incorporating  proper  bound¬ 
ary  conditions,  for  solving  the  Dirac  equation  in  arbitrary 
geometrical  confinements  in  (2  +  1)  dimensions.  Different 
from  the  k-p  models  (with  different  flavors  of  the  6-band 
or  8-band  Luttinger-Kohn  model),  a  simple  lattice  grid 
is  used  to  discretize  the  two-dimensional  closed  space. 
However,  in  solving  massless  Dirac  fermion  systems,  a 
necessary  ingredient  designed  to  eliminate  the  fermion 
doubling  problem  is  to  evaluate  the  Dirac  equations  and 
the  Dirac  spinors  at  two  different  sets  of  square  lattice 
points.  Specifically,  if  one  evaluates  the  Dirac  spinors  at 
locations  (z,  j),  the  Dirac  equations  have  to  be  evaluated  at 
(i  +  1/2,  j  +  1/2),  otherwise  one  ends  up  with  two  fermions 
in  the  first  Brillouin  zone,  which  conflicts  with  the  original 
single-fermion  Dirac  equation. 

For  non-relativistic  quantum  tunneling,  due  to  geomet¬ 
ric  symmetry  in  the  double- well  systems,  the  eigenstates 
are  either  symmetric  or  antisymmetric  about  the  central 
barrier.  In  this  case,  the  tunneling  rate  is  simply  propor¬ 
tional  to  the  energy  splitting  A E  between  the  pairing 
symmetric  and  antisymmetric  eigenstates  [4].  However, 
for  relativistic  quantum  tunneling,  only  in  one  spatial 
dimension  can  the  symmetric  and  antisymmetric  eigen¬ 
state  pairs  be  defined.  In  fact,  for  both  Dirac  fermion  in 
(2  +  1)  dimensions  and  graphene  systems,  the  symmet¬ 
ric  and  antisymmetric  eigenstates  are  not  necessarily  in 
pair.  A  more  general,  physically  meaningful  definition  of 
the  tunneling  rate  is  thus  needed.  Our  approach  is  the 
following.  For  an  arbitrary  symmetric  double-well  system, 
first  we  choose  a  random  linear  combination  of  eigenstates, 
denoted  as  0  =  anV+(^)-  We  then  set  values  of  0  on 

the  right  side  and  barrier  region  to  be  zero  and  renormal¬ 
ize  it.  The  new  state  is  denoted  as  0  =  an0^(r),  where 

an’s  are  the  renormalized  coefficients,  0^(r)  =  0n(r),  for 
r  in  the  left  well,  and  0^(r)  =0,  otherwise.  Next,  we  let 
this  state  0  evolve  with  time  and  express  it  in  terms  of  the 
linear  combination  of  all  eigenstates:  0  =  frm0m(^)5 

where  the  index  m  runs  through  all  eigenstates.  The 
coefficient  bm  can  be  calculated  as  bm  =  Id  V’mV’dr  = 
/l  E«  «++ =  En  /z,  V’mV’ndr,  where  D  and  L 
denote  the  integration  domains  of  the  whole  double  well 
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Fig.  2:  (Color  online)  Tunneling  rates  and  LDS  patterns  for 
massless  Dirac  fermion  in  integrable  and  chaotic  double- well 
systems,  where  0  and  %  are  two  components  of  the  Dirac  spinor. 
A  unit  system  h  =  c  =  1  has  been  used  in  our  calculations.  The 
barrier  height  Vo  is  about  60  under  such  unit  system.  For  the 
rectangle  double  well,  the  width  and  the  height  are  set  to 
2  and  1,  respectively.  The  bow-tie  chaotic  shape  is  obtained 
by  cutting  the  rectangle  with  three  arcs  such  that  the  cut 
parts  are  all  0.3  measured  on  both  sides  of  the  baseline  and 
central  vertical  line  of  the  rectangle.  The  theoretical  ratio  of 
the  left- well  width  to  the  barrier  width  should  be  24  :  1  for  the 
rectangle,  however,  this  ratio  may  vary  after  discretization. 


and  the  left  well,  respectively.  The  time  evolution  of 
the  state  ^  is  then  given  by  $(£)  =  brnil;rne~'lErnt^h. 
Because  the  particle  state  is  initially  confined  within  the 
left  well,  to  characterize  the  tunneling  process  of  this  state, 
we  calculate  the  probability  that  the  particle  may  be  found 
in  the  left  well  with  respect  to  time,  P/,(t)  =  fL  \ip(t)\2dr. 
The  tunneling  rate  R  can  be  determined  when  Pl^) 
reaches  the  minimum  for  the  first  time,  i.e.,R  =  7rAP/AT, 
where  A P  is  the  probability  difference  between  the  initial 
value  and  the  first  minimum  of  Pl(£),  and  AT  is  the 
time  it  takes  to  reach  the  minimum.  This  definition  is 
general  because  for  non-relativistic  quantum  tunneling 
where  symmetric/antisymmetric  eigenstate  pairs  do  exist, 
it  reduces  to  A E. 

We  now  present  numerical  evidence  for  the  effect  of 
chaos  on  relativistic  quantum  tunneling.  Figure  2  shows 
the  generalized  tunneling  rate  R  vs.  the  normalized  energy 
E/Vo  for  massless  Dirac  fermion  in  the  double- well  barrier 
system  for  two  types  of  geometry:  one  classically  inte¬ 
grable  and  another  chaotic.  For  the  integrable  geometry, 
we  observe  the  existence  of  states  with  extremely  low 
tunneling  rates,  as  indicated  by  the  arrow  in  fig.  2(a). 
These  correspond  to  states  localized  nearly  entirely  in  the 
left  or  right  side  of  the  potential  barrier,  which  “survive” 
the  tunneling  process  between  the  two  sides,  as  indicated 
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Fig.  3:  (Color  online)  Tunneling  rates  and  LDS  patterns  in 
integrable  and  chaotic  graphene  double- well  systems,  where  A 
and  B  denote  the  two  distinct  types  of  atoms  in  a  graphene 
unit  cell.  The  tight-binding  model  is  used  in  the  calculation 
for  graphene  systems.  The  barrier  height  is  fixed  at  Vo  ~  0.67t 
(with  t  denoting  the  nearest-neighbor  hopping  integral).  After 
converting  to  the  h  =  c  =  1  unit  system,  the  geometric  measures 
and  the  barrier  potential  are  the  same  as  the  single  Dirac 
fermion  system. 


by  the  accompanying  pattern  of  local  density  of  states 
(LDS).  We  note  that,  in  non-relativistic  quantum  trans¬ 
port,  these  are  effectively  quantum  pointer  states  [13,14]. 
In  relativistic  quantum  systems,  we  observe  that  both 
components  of  the  underlying  Dirac  spinor  exhibit  a  heavy 
concentration  of  the  probability  in  orbits  along  which 
particles  travel  vertically  back  and  forth  on  either  side, 
parallel  to  the  barrier.  For  the  chaotic  geometry,  while 
signatures  of  pointer  states  can  still  be  found,  they  are 
weak  as  compared  with  those  in  the  integrable  counter¬ 
part,  as  shown  in  fig.  2(b)  and  the  accompanying  LDS 
pattern.  Analogous  behaviors  occur  when  the  entire  cavity 
is  made  of  graphene,  as  shown  in  fig.  3.  Thus,  in  both 
Dirac  fermion  and  graphene  systems,  we  observe  that  clas¬ 
sical  chaos  can  greatly  suppress  the  spread  in  the  quantum 
tunneling  rate,  as  in  non-relativistic  quantum  systems  [4]. 

A  common  phenomenon  between  Dirac-fermion  and 
graphene  tunneling  systems  is  that,  for  small  energies, 
pointer  states  are  far  less  likely  than  in  non-relativistic 
quantum  systems.  Consequently,  in  both  systems,  the 
tunneling  rate  can  be  quite  large  even  in  the  small  energy 
regime,  as  shown  in  figs.  2  and  3.  This  is  the  direct 
consequence  of  Klein  tunneling  [6-8],  which  finds  no 
counterpart  in  non-relativistic  quantum  mechanics  where 
the  tunneling  rate  tends  to  zero  as  the  energy  is  decreased 
to  zero  [4].  Although  both  systems  show  almost  identical 
relativistic  behaviors  at  low  energies,  and  both  prove 
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chaos  regularization  to  be  a  universal  phenomenon  across 
quantum  systems,  one  should  not  confuse  one  system  with 
the  other.  The  difference  between  the  two  systems  lies  in 
the  number  of  massless  (quasi-)particles.  It  is  clear  that 
the  system  described  by  eq.  (1)  contains  a  single  massless 
Dirac  fermion.  On  the  other  hand,  for  a  graphene  system 
in  the  low-energy  limit,  the  two  electron  states  at  different 
atoms  (commonly  called  A  and  B)  in  a  unit  cell  can  be 
thought  of  as  the  two  states  of  a  massless  quasi-particle’s 
pseudo-spin.  However,  there  are  actually  two  of  these 
massless  particles  close  to  the  two  distinct  Dirac  points. 
In  the  presence  of  a  spatially  short-range  potential,  the 
coupling  between  the  two  Dirac  points  becomes  apparent. 

To  understand  the  effect  of  chaos  on  relativistic  quan¬ 
tum  tunneling,  we  develop  a  theory  based  on  the  self¬ 
energy  concept  widely  employed  in  the  study  of  quantum 
transport  [15,16].  The  basic  observation  is  that  pointer 
states  generally  result  in  low  tunneling  rate.  Thus,  when 
pointer  states  are  present,  the  coupling  between  the  two 
wells  is  weak.  The  tunneling  rate  can  thus  be  approxi¬ 
mated  as  the  escaping  rate  of  Dirac  fermions  between  two 
nearly  closed,  weakly  coupled  wells.  Let  the  left  well  be 
denoted  by  superscript  (1),  and  the  barrier  together  with 
the  right  well  denoted  by  superscript  (2).  The  Dirac  equa¬ 
tions  for  the  whole  double-well  system  can  be  written  in 
terms  of  the  Hamiltonians  Hi  and  H 2  for  the  separated 
closed  wells  as 


~HX  V12 

—  7 7 

xp  C1) 

V21  H-2 

*0(2) 

—  Hz 

-0(2) 

where  Vij  are  the  coupling  matrices.  Note  that,  if 
the  left  well  were  itself  closed,  the  corresponding 
equation  would  be  HiipW  =  Eixp^\  which  becomes 
(Hi  +  =  Eip^  when  weak  coupling  to  the  right 

well  is  taken  into  account,  where  T,R  =  V\2GRV2 1  is  the 
self-energy  due  to  the  barrier  and  the  right  well,  and 
GR  =  (E  +  ir1-H2)-1  is  the  retarded  Green’s  function. 
For  each  eigenstate  in  the  left  well,  the  energy  shift  can 
be  obtained  through  first-order  perturbation  theory  as 
(E^)  =  (ip^ |SjR|/0(1)),  which  is  typically  complex.  The 
real  part  of  this  shift  changes  the  oscillating  frequency  of 
the  corresponding  eigenstate,  while  the  imaginary  part, 
denoted  by  —7,  introduces  a  decay  factor  exp(— cojt/h) 
in  the  time  evolution  of  the  probability,  which  describes 
the  escaping  rate  of  the  Dirac  fermions  from  the  left 
to  the  right  well.  Note  that,  since  the  whole  system  is 
closed,  7  only  describes  the  transient  behavior  associated 
with  particle’s  tunneling  from  left  to  right,  with  any 
recurring  behavior  neglected.  If  we  let  the  right  well 
extend  to  infinity  so  that  there  is  no  reflection,  the 
system  is  equivalent  to  a  single  left  well  coupled  with  a 
semi-infinite  lead  through  a  potential  barrier,  and  7  will 
be  the  tunneling  rate  for  the  single-left-well  system. 

To  gain  insights  and  also  to  validate  that  our  proposed 
quantity  7  is  indeed  proportional  to  the  tunneling  rate, 
we  consider  a  one-dimensional  system  for  which  the 
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Fig.  4:  (Color  online)  Schematic  diagram  for  one-dimension 
tunneling. 


Dirac  equation  can  be  solved  and  the  self-energy  and 
consequently  7  can  then  be  calculated  analytically  through 
the  Green’s  function,  as  shown  in  fig.  4.  Specifically,  the 
one-dimensional  Dirac  equation  is  (—ihvaxdx)xp  =  (E  — 
V)ip.  Since  the  potential  is  zero  in  the  left  well  (of  width 
Lx),  the  solution  is 


tf)(x>  =  7Cexp(il) 

where  kn  =  (n-\-l/2)ir/Lx.  For  the  right  part,  the  barrier 
has  width  W  and  V  =  Vo ,  and  the  right  well  has  width  L 
and  V  =  0.  We  obtain  the  solution 


cos (knx-  f ) 
ism(knx-j) 


V> {n\x ) 


AielKnXu+  -\-A2e  1KnXU-,  barrier, 
A^elknXu+  +  A^e~lknXU- ,  right  well, 


where  kn  =  E/(hv)  =  [V0W/(hv)  +  (m  +  1/2)tt]/(L  +  W), 
Kn  =  kn  —  Vb/(^u),  u+  =  (1, 1)T  and  u- =  ( 1,  — 1)T  are 
bases  for  the  spinor,  and  the  coefficients  Ai  are  determined 
by  the  boundary  conditions.  Note  that  this  is  a  combined 
solution  for  both  E  <  Vo  and  E>Vq.  To  calculate  the 
self-energy  (E^),  we  discretize  the  space  and  introduce 
the  Green’s  function  GR(xi,Xj).  The  self-energy  can  be 
expressed  as  [15] 

(Sfi)  =  ^(Lx-a)V12GR(Lx  +a,Lx  +  a)V21^(Lx-a), 


where  x  =  Lx  is  the  junction  between  the  left  well  and 
the  barrier,  and  the  coupling  components  are  V\2  =  = 

—ihvax.  The  Green’s  function  can  be  calculated  explicitly 
in  the  limit  L  — >•  00.  After  a  substantial  amount  of  algebra 
we  obtain  the  following  expression  for  the  tunneling  rate: 

7  =  -  [cos(2fca)  Im(C')  -  sin(2A;a)  In^S)]  -  , 

where  the  quantities  S  and  C  are  given  by 

sin 

E  +  iri  —  hvk 5 
cos(2  KCl 
E  +  iri  ~ 
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Fig.  5:  (Color  online)  Theoretical  tunneling  rate  7  for  Dirac 
fermion  in  the  two-dimensional  ring  system  as  described  in 
the  text.  The  eigenstates  associated  with  the  inner  and  outer 
ring  regions  are  calculated  numerically  by  solving  the  Dirac 
equation  with  proper  boundary  conditions.  The  parameters  are 
Ri  =  5,  R2  =  10,  and  Vo  =  5.  The  dimensions  for  the  outer  part 
are  R3  =  10.2,  and  R4  —  R3  =  3(i?2  —  Ri). 

The  integrals  can  be  evaluated  explicitly  and  we  finally 
obtain  y  &hv/(2Lx).  In  one  spatial  dimension  the  stan¬ 
dard  tunneling  rate  A E  can  also  be  analytically  calcu¬ 
lated  through  symmetry  considerations.  We  obtain  A E  ~ 
Hv7t/(2Lx)  and  consequently,  7  «  AE / 7r.  These  analytic 
predictions  have  been  verified  numerically. 

For  tunneling  systems  in  two  spatial  dimensions,  even 
when  the  geometry  is  regular  it  is  not  possible  to  calcu¬ 
late  the  tunneling  rate  analytically,  due  to  the  entan¬ 
glement  of  the  two  Cartesian  coordinates  in  the  Dirac 
equation.  The  problem  is  solvable  only  for  certain  special 
types  of  boundary  conditions  for  which  the  variables  can 
be  separated.  One  particular  class  of  such  systems  is 
tunneling  between  two  concentric  ring  regions,  where  the 
first  cavity  is  given  by  Ri  <  r  <  R2,  the  potential  barrier 
is  located  in  the  region  R2  <r  <  R%,  and  the  second  cavity 
is  in  the  region  R3  <  r  <  R4.  When  all  four  radii  are  large, 
this  ring-shaped  double- well  tunneling  system  is  topologi¬ 
cally  equivalent  to  a  rectangular  double  well  with  periodic 
boundary  condition  in  one  direction.  We  have  solved  the 
Dirac  equation,  calculated  the  Green’s  function  and  the 
self-energy,  and  finally  obtained  a  close-form  expression 
for  the  tunneling  rate  7,  which  involves  various  eigenstates 
that  can  be  evaluated  numerically.  A  representative  behav¬ 
ior  of  7  as  a  function  of  the  energy  is  shown  in  fig.  5,  which 
is  qualitatively  similar  to  that  in  the  rectangular  double¬ 
well  system.  Of  particular  importance  to  the  aim  of  this 
letter  are  pointer  states  with  extraordinarily  low  tunneling 
rates,  which  are  the  causes  of  the  tunneling  spread. 

In  summary,  we  have  established  the  phenomenon  of 
chaos-regularized  tunneling  in  the  relativistic  quantum 
regime  for  both  Dirac  fermion  and  graphene  tunnel¬ 
ing  systems.  While  tunneling  in  non-relativistic  quantum 


mechanics  has  been  well  studied,  there  are  two  major 
challenges  in  investigating  the  phenomenon  in  relativistic 
quantum  systems:  solution  of  the  Dirac  equation  in  general 
two-dimensional  spatial  geometry  and  physically  mean¬ 
ingful  definition  of  the  tunneling  rate.  We  have  overcome 
these  difficulties  and  developed  numerical  methods  and 
analytic  theory  to  reveal  the  generic  features  of  relativistic 
quantum  tunneling.  While  we  focus  on  resonant  tunneling, 
the  methodology  developed  here  can  be  extended  to  other 
relativistic  quantum  transport  systems,  such  as  various 
electronic  transport  devices  made  of  graphene. 
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